/**
 * @file   xbandec.c
 * @author  <xujian@xj-workstation>
 * @date   Thu Jan 12 16:44:52 2012
 * 
 * @brief  It is an example from NR
 *         float -> double replaced by me
 * 
 *         Solving the 1D-NLD with v != 0 
 *         split-step scheme
 *         3-order split error ?
 */


/* Driver for routine bandec */

#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <time.h>
#include "nr.h"
#include "nrutil.h"
#include <pthread.h>

#define NRANSI

int splitstep1(double *x, double *xs, int Nx,
			   double dx, double pa_m, double dt,
			   int thread_number, int flag);
int splitstep2(double *x, double *xs, int Nx,
			   double pa_m, double pa_lam,
			   double dt, int thread_number, int flag);
/** 
 * 
 * 
 * @param x 
 * @param xs 
 * @param Nx 
 * @param dx 
 * @param pa_m 
 * @param dt 
 * @param thread_number 
 * @param n_step 
 * 
 * @return xs
 */
int splitFoward(double *x, double *xs, int Nx, double dx,
				double pa_m, double pa_lam,double dt,
				int thread_number, int n_step){
		splitstep1(x, xs, Nx, dx, pa_m,
				   dt, thread_number, n_step);		
		splitstep2(xs, x, Nx, pa_m, pa_lam, dt,
				   thread_number, n_step);
		return 0;
}
/** 
 * 
 * 
 * @param x 
 * @param xs 
 * @param Nx 
 * @param dx 
 * @param pa_m 
 * @param dt 
 * @param thread_number 
 * @param n_step 
 * 
 * @return xs
 */
int  splitadjoint(double *x, double *xs, int Nx, double dx,
				  double pa_m, double pa_lam,double dt,
				  int thread_number, int n_step){
		splitstep2(x, xs, Nx, pa_m, pa_lam, dt,
				   thread_number, n_step);
		splitstep1(xs, x, Nx, dx, pa_m,
				   dt, thread_number, n_step);
		return 0;
}





double dcharge(double *x, int Nx, double dx);
double denergy(double *x, int Nx, double dx, double pa_m, double pa_lam);
double dmomentum(double *x, int Nx, double dx);
int accurary(double *x, double *rx, int Nx, double dx,
			 double *L_inf, double *R_inf);
int output(double q_phy[3][3], double *x, double *rx, double *rtx,
		   int Nx, double dx, double pa_m, double pa_lam, double pa_La[],
		   double dt, double time, int n_cyc, FILE *fp, double xlr[], double Xl);
int initialization(double **a, double **al, double **ar, double q_phy[3][3],
				   double *x, double *rx, int Nx, double dx, double pa_m,
				   double pa_lam, double lam4, double pa_La[], double Xl,
				   double Tmax, double xlr[]);
/* int seta(double **a, double **ar, double *x, double *xs, int Nx, */
/* 		 double lam4, double pa_m, double pa_lam, double dt); */
int	outputu(double *x, double *rx, int Nx, int n_step, double dt, FILE **fpu);
double loneerr1(double *a, double *b, int Nx);
double loneerr2(double *a, double *b, int Nx);
int ABx(double pa_m, double pa_lam, double pa_La, double x,
		double *tempA, double *tempB);
int
copya2b(double *a, double *b, int Nr){
	int i;
	for(i = 1; i <= Nr; i++)
		b[i] = a[i];
	return 0;
}

int
sign(double x){
	if(x>1e-14)
		return 1;
	else if(x<1e-14)
		return -1;
	else
		return 0;
}			

int main(void)
{
	FILE *fp;					/**< 文件输入指针 */
	char s[255];				/**< 文件读取注释 */
	FILE **fpu, *fpe;			/**< 文件输出指针 */
    int err;					/**< 函数错误返回值 */
	double ep = 1e-10;			/**< 迭代允许误差 */
    double d,**a,**al,*b,*x;    /**< 求解带状阵所需变量 */
	double **ar;				/**< 等式右端带状矩阵 */
    double *rtx, *rx, *xs;	    /**< 真解向量, 真解未乘时间
								   \f$ exp{-i\Lambda t} \f$,迭代时储存向量 */
    int Nx, Nt, Nr;	            /**< 空间网格数,时间步数,储存向量长度 */
    double Xl, Xr;		        /**< 空间左右界限 */
	double xlr[2];				/**< \f$ x_0 \f$ and velocity*/
    double Tmax;		        /**< 时间终点 */
    double dt, dx, lambda;		/**< 时间,空间网格步长,dt/dx */
    unsigned long *indx;        /**< 排列阵指针*/
	int  i, n_step, n_cyc;	    /**< 循环变量 */
    double pa_m, pa_La[2], pa_lam;	/**< parameter */
    double lam4;
	/* double *tempp; */
    double q_phy[3][3];	/**< [charge; energy; momentum]
						 * [初值，数值解，相对误差] */
	int maxitera, minitera;		/**< 最大，最小迭代步数 */
	int outputatstep;			/**< 每隔多少步输出 */
	double time;
	
	clock_t start;
	clock_t end;
	start = clock();
	
	int thread_number;

	/* pthread_t *thr_id; */
	/* struct arg_t *args; */
	

	/* pthread_mutex_t lock; */
	int as[18] = {1,1,1,-2,1,1,1,1,1};
	int a_flag[18] = {0,1,0,1,0,0,0,0,1};
	for(i=0; i<9; i++){
		as[17-i] = as[i];
		a_flag[17-i] = !a_flag[i];
	}
	
	
	fp = fopen("parameter", "r");
	fpu = (FILE **) malloc(2*sizeof(FILE*));
	fpu[0] = fopen("numericalS", "w");
	fpu[1] = fopen("iteration", "w");
	fpe = fopen("error", "w");
	if(!fp){
		printf("fail to open the file <parameter>\n");
		return 1;
	}
	else{
		if(fscanf(fp,"%lf", &Tmax))
			fprintf(fpu[1],"%f", Tmax);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%d", &Nx))
			fprintf(fpu[1],"%d", Nx);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%d", &Nt))
			fprintf(fpu[1],"%d", Nt);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%lf", &ep))
			fprintf(fpu[1],"%e", ep);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		
		if(fscanf(fp,"%lf", &pa_m))
			fprintf(fpu[1],"%f", pa_m);
		else return 1;
		fprintf(fpu[1],"%s", s);
		fgets(s, 255, fp);
		if(fscanf(fp,"%lf", &pa_lam))
			fprintf(fpu[1],"%f", pa_lam);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);

		if(fscanf(fp,"%lf", xlr))
			fprintf(fpu[1],"%f", xlr[0]);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%lf", xlr+1))
			fprintf(fpu[1],"%f", xlr[1]);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);

		if(fscanf(fp,"%lf", pa_La))
			fprintf(fpu[1],"%f", pa_La[0]);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%lf", pa_La+1))
			fprintf(fpu[1],"%f", pa_La[1]);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		
		if(fscanf(fp,"%lf", &Xl))
			fprintf(fpu[1],"%f", Xl);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%lf", &Xr))
			fprintf(fpu[1],"%f", Xr);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);		
		if(fscanf(fp,"%d", &maxitera))
			fprintf(fpu[1],"%d", maxitera);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
		if(fscanf(fp,"%d", &minitera))
			fprintf(fpu[1],"%d", minitera);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);

		if(fscanf(fp,"%d", &outputatstep ))
			fprintf(fpu[1],"%d", outputatstep);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);

		if(fscanf(fp,"%d", &thread_number ))
			fprintf(fpu[1],"%d", thread_number);
		else return 1;
		fgets(s, 255, fp);
		fprintf(fpu[1],"%s", s);
	}
		
    dx = (Xr - Xl)/(Nx-1);	
    Nr = 4*Nx;	
    dt = Tmax/Nt;
	lambda = dt/dx;

	lam4 = lambda/4.0;

	
    /// allocte space to the matrix and pointers
    a = dmatrix(1,Nr,1,13);
	ar = dmatrix(1,Nr,1,13);
    x = dvector(1,Nr);
    b = dvector(1,Nr);
    al = dmatrix(1,Nr,1,6);
    indx = lvector(1,Nr);
    rx = dvector(1,Nr);
    xs = dvector(1,Nr);
    rtx = dvector(1,Nr);
	err = initialization(a, al, ar, q_phy, x, rtx, Nx, dx, pa_m,
						 pa_lam, lam4, pa_La, Xl, Tmax, xlr);
	
	err = copya2b(x, xs, Nr);

	/* fprintf(fpu[0], "%16.5f  ", dt); */
	/* for(j=1;j<=4;j++){ */
	/* 	for(i=0; i<Nx; i++){ */
	/* 		fprintf(fpu[0], "%18.14e  ", Xl+i*dx); */
	/* 	} */
	/* } */
	/* fprintf(fpu[0], "\n"); */

	n_step = 0;
	n_cyc = 0;
	time = 0.0;
	while(time < Tmax-dt){
		/* splitadjoint(x, xs, Nx, dx, pa_m, pa_lam, dt, */
		/* 			thread_number, n_step); */
		/* for(i=0; i<3; i++){ */
		/* 			splitFoward(x, xs, Nx, dx, pa_m, pa_lam, dt, */
		/* 						thread_number, n_step); */
		/* } */
		/* splitadjoint(x, xs, Nx, dx, pa_m, pa_lam, dt, */
		/* 			thread_number, n_step); */
		/* splitadjoint(x, xs, Nx, dx, pa_m, pa_lam, -2*dt, */
		/* 			thread_number, n_step); */
		/* for(i=0; i<3; i++){ */
		/* 			splitFoward(x, xs, Nx, dx, pa_m, pa_lam, dt, */
		/* 						thread_number, n_step); */
		/* }				 */

		for(i=0; i<18; i++){
			if(a_flag[i]){
				splitFoward(x, xs, Nx, dx, pa_m, pa_lam, as[i]*dt,
							thread_number, n_step);				
			}
			else{
				splitadjoint(x, xs, Nx, dx, pa_m, pa_lam, as[i]*dt,
							 thread_number, n_step);
			}
			time += as[i]*dt;
		}
		/* time += 6*dt; */
		output(q_phy, xs, rx, rtx, Nx, dx, pa_m, pa_lam,
			   pa_La, dt, time, n_cyc, fpe, xlr, Xl);

	}
	output(q_phy, xs, rx, rtx, Nx, dx, pa_m, pa_lam,
		   pa_La, dt, time, n_cyc, fpe, xlr, Xl);


/// strang split below
	/* splitstep2(xs, x, Nx, pa_m, pa_lam, dt/2, */
	/* 		   thread_number, n_step); */
	/* time = dt/2; */
	/* while(time < Tmax-dt){		 */
	/* 	splitstep1(x, xs, Nx, dx, pa_m, */
	/* 			   dt, thread_number, n_step); */
	/* 	splitstep2(xs, x, Nx, pa_m, pa_lam, dt, */
	/* 			   thread_number, n_step); */

	/* 	/\* tempp = xs; *\/ */
	/* 	/\* xs = x; *\/ */
	/* 	/\* x = tempp; *\/ */
	/* 	time += dt; */
	/* } */
	/* splitstep1(x, xs, Nx, dx, pa_m, */
	/* 		   dt, thread_number, n_step); */
	/* splitstep2(xs, x, Nx, pa_m, pa_lam, dt/2, */
	/* 		   thread_number, n_step); */
	/* output(q_phy, xs, rx, rtx, Nx, dx, pa_m, pa_lam, */
	/* 	   pa_La, dt, time+dt/2, n_cyc, fpe, xlr, Xl); */
	
    /// free the pointers
    free_dvector(b,1,Nr);
    free_dvector(x,1,Nr);
    free_dvector(rx,1,Nr);
    free_dvector(rtx,1,Nr);
    free_dvector(xs,1,Nr);
	fclose(fp);

	end = clock();
	fprintf(fpu[1],"\n cost time:%ld\n", (end - start)*1000/CLOCKS_PER_SEC);
	fclose(fpu[0]);
	fclose(fpu[1]);
	fclose(fpe);
	free(fpu);

	
	return 0;
}

#undef NRANSI
